home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMAT4.CPP < prev    next >
C/C++ Source or Header  |  1995-01-17  |  19KB  |  689 lines
  1. //$$ newmat4.cpp       Constructors, ReDimension, basic utilities
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10. //#define REPORT { static ExeCounter ExeCount(__LINE__,4); ++ExeCount; }
  11.  
  12. #define REPORT {}
  13.  
  14.  
  15.  
  16.  
  17.  
  18.  
  19.  
  20.  
  21. /*************************** general utilities *************************/
  22.  
  23. static int tristore(int n)                      // els in triangular matrix
  24. { return (n*(n+1))/2; }
  25.  
  26.  
  27. /****************************** constructors ***************************/
  28.  
  29. GeneralMatrix::GeneralMatrix()
  30. { store=0; storage=0; nrows=0; ncols=0; tag=-1; }
  31.  
  32. GeneralMatrix::GeneralMatrix(ArrayLengthSpecifier s)
  33. {
  34.    REPORT
  35.    storage=s.Value(); tag=-1;
  36.    if (storage)
  37.    {
  38.       store = new Real [storage]; MatrixErrorNoSpace(store);
  39.       MONITOR_REAL_NEW("Make (GenMatrix)",storage,store)
  40.    }
  41.    else store = 0;
  42. }
  43.  
  44. Matrix::Matrix(int m, int n) : GeneralMatrix(m*n)
  45. { REPORT nrows=m; ncols=n; }
  46.  
  47. SymmetricMatrix::SymmetricMatrix(ArrayLengthSpecifier n)
  48.    : GeneralMatrix(tristore(n.Value()))
  49. { REPORT nrows=n.Value(); ncols=n.Value(); }
  50.  
  51. UpperTriangularMatrix::UpperTriangularMatrix(ArrayLengthSpecifier n)
  52.    : GeneralMatrix(tristore(n.Value()))
  53. { REPORT nrows=n.Value(); ncols=n.Value(); }
  54.  
  55. LowerTriangularMatrix::LowerTriangularMatrix(ArrayLengthSpecifier n)
  56.    : GeneralMatrix(tristore(n.Value()))
  57. { REPORT nrows=n.Value(); ncols=n.Value(); }
  58.  
  59. DiagonalMatrix::DiagonalMatrix(ArrayLengthSpecifier m) : GeneralMatrix(m)
  60. { REPORT nrows=m.Value(); ncols=m.Value(); }
  61.  
  62. Matrix::Matrix(const BaseMatrix& M)
  63. {
  64.    REPORT // CheckConversion(M);
  65.    MatrixConversionCheck mcc;
  66.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Rt);
  67.    GetMatrix(gmx);
  68. }
  69.  
  70. RowVector::RowVector(const BaseMatrix& M) : Matrix(M)
  71. {
  72.    if (nrows!=1)
  73.    {
  74.       Tracer tr("RowVector");
  75.       Throw(VectorException(*this));
  76.    }
  77. }
  78.  
  79. ColumnVector::ColumnVector(const BaseMatrix& M) : Matrix(M)
  80. {
  81.    if (ncols!=1)
  82.    {
  83.       Tracer tr("ColumnVector");
  84.       Throw(VectorException(*this));
  85.    }
  86. }
  87.  
  88. SymmetricMatrix::SymmetricMatrix(const BaseMatrix& M)
  89. {
  90.    REPORT  // CheckConversion(M);
  91.    MatrixConversionCheck mcc;
  92.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Sm);
  93.    GetMatrix(gmx);
  94. }
  95.  
  96. UpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
  97. {
  98.    REPORT // CheckConversion(M);
  99.    MatrixConversionCheck mcc;
  100.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UT);
  101.    GetMatrix(gmx);
  102. }
  103.  
  104. LowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
  105. {
  106.    REPORT // CheckConversion(M);
  107.    MatrixConversionCheck mcc;
  108.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LT);
  109.    GetMatrix(gmx);
  110. }
  111.  
  112. DiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
  113. {
  114.    REPORT //CheckConversion(M);
  115.    MatrixConversionCheck mcc;
  116.    GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Dg);
  117.    GetMatrix(gmx);
  118. }
  119.  
  120. GeneralMatrix::~GeneralMatrix()
  121. {
  122.    if (store)
  123.    {
  124.       MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
  125. #ifdef Version21
  126.       delete [] store;
  127. #else
  128.       delete [storage] store;
  129. #endif
  130.    }
  131. }
  132.  
  133. CroutMatrix::CroutMatrix(const BaseMatrix& m)
  134. {
  135.    REPORT
  136.    Tracer tr("CroutMatrix");
  137.    GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::Rt);
  138.    GetMatrix(gm);
  139.    if (nrows!=ncols) Throw(NotSquareException(*this));
  140.    d=TRUE; sing=FALSE;
  141.    indx=new int [nrows]; MatrixErrorNoSpace(indx);
  142.    MONITOR_INT_NEW("Index (CroutMat)",nrows,indx)
  143.    ludcmp();
  144. }
  145.  
  146. CroutMatrix::~CroutMatrix()
  147. {
  148.    MONITOR_INT_DELETE("Index (CroutMat)",nrows,indx)
  149. #ifdef Version21
  150.    delete [] indx;
  151. #else
  152.    delete [nrows] indx;
  153. #endif
  154. }
  155.  
  156. //ReturnMatrixX::ReturnMatrixX(GeneralMatrix& gmx)
  157. //{
  158. //   REPORT
  159. //   gm = gmx.Image(); gm->ReleaseAndDelete();
  160. //}
  161.  
  162. #ifndef TEMPS_DESTROYED_QUICKLY_R
  163.  
  164. GeneralMatrix::operator ReturnMatrixX() const
  165. {
  166.    REPORT
  167.    GeneralMatrix* gm = Image(); gm->ReleaseAndDelete(); 
  168.    return ReturnMatrixX(gm);
  169. }
  170.  
  171. #else
  172.  
  173. GeneralMatrix::operator ReturnMatrixX&() const
  174. {
  175.    REPORT
  176.    GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
  177.    ReturnMatrixX* x = new ReturnMatrixX(gm);
  178.    MatrixErrorNoSpace(x); return *x;
  179. }
  180.  
  181. #endif
  182.  
  183. #ifndef TEMPS_DESTROYED_QUICKLY_R
  184.  
  185. ReturnMatrixX GeneralMatrix::ForReturn() const
  186. {
  187.    REPORT
  188.    GeneralMatrix* gm = Image(); gm->ReleaseAndDelete(); 
  189.    return ReturnMatrixX(gm);
  190. }
  191.  
  192. #else
  193.  
  194. ReturnMatrixX& GeneralMatrix::ForReturn() const
  195. {
  196.    REPORT
  197.    GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
  198.    ReturnMatrixX* x = new ReturnMatrixX(gm);
  199.    MatrixErrorNoSpace(x); return *x;
  200. }
  201.  
  202. #endif
  203.  
  204. /**************************** ReDimension matrices ***************************/
  205.  
  206. void GeneralMatrix::ReDimension(int nr, int nc, int s)
  207. {
  208.    REPORT 
  209.    if (store)
  210.    {
  211.       MONITOR_REAL_DELETE("Free (ReDimensi)",storage,store)
  212. #ifdef Version21
  213.       delete [] store;
  214. #else
  215.       delete [storage] store;
  216. #endif
  217.    }
  218.    storage=s; nrows=nr; ncols=nc; tag=-1;
  219.    if (s)
  220.    {
  221.       store = new Real [storage]; MatrixErrorNoSpace(store);
  222.       MONITOR_REAL_NEW("Make (ReDimensi)",storage,store)
  223.    }
  224.    else store = 0;
  225. }
  226.  
  227. void Matrix::ReDimension(int nr, int nc)
  228. { REPORT GeneralMatrix::ReDimension(nr,nc,nr*nc); }
  229.  
  230. void SymmetricMatrix::ReDimension(int nr)
  231. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  232.  
  233. void UpperTriangularMatrix::ReDimension(int nr)
  234. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  235.  
  236. void LowerTriangularMatrix::ReDimension(int nr)
  237. { REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  238.  
  239. void DiagonalMatrix::ReDimension(int nr)
  240. { REPORT GeneralMatrix::ReDimension(nr,nr,nr); }
  241.  
  242. void RowVector::ReDimension(int nc)
  243. { REPORT GeneralMatrix::ReDimension(1,nc,nc); }
  244.  
  245. void ColumnVector::ReDimension(int nr)
  246. { REPORT GeneralMatrix::ReDimension(nr,1,nr); }
  247.  
  248. void RowVector::ReDimension(int nr, int nc)
  249. {
  250.    Tracer tr("RowVector::ReDimension");
  251.    if (nr != 1) Throw(VectorException(*this));
  252.    REPORT GeneralMatrix::ReDimension(1,nc,nc);
  253. }
  254.  
  255. void ColumnVector::ReDimension(int nr, int nc)
  256. {
  257.    Tracer tr("ColumnVector::ReDimension");
  258.    if (nc != 1) Throw(VectorException(*this));
  259.    REPORT GeneralMatrix::ReDimension(nr,1,nr);
  260. }
  261.  
  262.  
  263. /********************* manipulate types, storage **************************/
  264.  
  265. int GeneralMatrix::search(const BaseMatrix* s) const
  266. { REPORT return (s==this) ? 1 : 0; }
  267.  
  268. int GenericMatrix::search(const BaseMatrix* s) const
  269. { REPORT return gm->search(s); }
  270.  
  271. int MultipliedMatrix::search(const BaseMatrix* s) const
  272. { REPORT return bm1->search(s) + bm2->search(s); }
  273.  
  274. int ShiftedMatrix::search(const BaseMatrix* s) const
  275. { REPORT return bm->search(s); }
  276.  
  277. int NegatedMatrix::search(const BaseMatrix* s) const
  278. { REPORT return bm->search(s); }
  279.  
  280. int ConstMatrix::search(const BaseMatrix* s) const
  281. { REPORT return (s==cgm) ? 1 : 0; }
  282.  
  283. int ReturnMatrixX::search(const BaseMatrix* s) const
  284. { REPORT return (s==gm) ? 1 : 0; }
  285.  
  286. MatrixType Matrix::Type() const { return MatrixType::Rt; }
  287. MatrixType SymmetricMatrix::Type() const { return MatrixType::Sm; }
  288. MatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }
  289. MatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }
  290. MatrixType DiagonalMatrix::Type() const { return MatrixType::Dg; }
  291. MatrixType RowVector::Type() const { return MatrixType::RV; }
  292. MatrixType ColumnVector::Type() const { return MatrixType::CV; }
  293. MatrixType CroutMatrix::Type() const { return MatrixType::Ct; }
  294. MatrixType BandMatrix::Type() const { return MatrixType::BM; }
  295. MatrixType UpperBandMatrix::Type() const { return MatrixType::UB; }
  296. MatrixType LowerBandMatrix::Type() const { return MatrixType::LB; }
  297. MatrixType SymmetricBandMatrix::Type() const { return MatrixType::SB; }
  298.  
  299. MatrixBandWidth BaseMatrix::BandWidth() const { return -1; }
  300. MatrixBandWidth DiagonalMatrix::BandWidth() const { return 0; }
  301.  
  302. MatrixBandWidth BandMatrix::BandWidth() const
  303.    { return MatrixBandWidth(lower,upper); }
  304.  
  305. MatrixBandWidth GenericMatrix::BandWidth() const { return gm->BandWidth(); }
  306.  
  307. MatrixBandWidth AddedMatrix::BandWidth() const
  308. { return gm1->BandWidth() + gm2->BandWidth(); }
  309.  
  310. MatrixBandWidth SPMatrix::BandWidth() const
  311. { return gm1->BandWidth().minimum(gm2->BandWidth()); }
  312.  
  313. MatrixBandWidth MultipliedMatrix::BandWidth() const
  314. { return gm1->BandWidth() * gm2->BandWidth(); }
  315.  
  316. MatrixBandWidth ConcatenatedMatrix::BandWidth() const { return -1; }
  317. MatrixBandWidth SolvedMatrix::BandWidth() const { return -1; }
  318. MatrixBandWidth ScaledMatrix::BandWidth() const { return gm->BandWidth(); }
  319. MatrixBandWidth NegatedMatrix::BandWidth() const { return gm->BandWidth(); }
  320.  
  321. MatrixBandWidth TransposedMatrix::BandWidth() const
  322. { return gm->BandWidth().t(); }
  323.  
  324. MatrixBandWidth InvertedMatrix::BandWidth() const { return -1; }
  325. MatrixBandWidth RowedMatrix::BandWidth() const { return -1; }
  326. MatrixBandWidth ColedMatrix::BandWidth() const { return -1; }
  327. MatrixBandWidth DiagedMatrix::BandWidth() const { return 0; }
  328. MatrixBandWidth MatedMatrix::BandWidth() const { return -1; }
  329. MatrixBandWidth ConstMatrix::BandWidth() const { return cgm->BandWidth(); }
  330. MatrixBandWidth ReturnMatrixX::BandWidth() const { return gm->BandWidth(); }
  331.  
  332. MatrixBandWidth GetSubMatrix::BandWidth() const
  333. {
  334.  
  335.    if (row_skip==col_skip && row_number==col_number) return gm->BandWidth();
  336.    else return MatrixBandWidth(-1);
  337. }
  338.  
  339. /************************ the memory managment tools **********************/
  340.  
  341. //  Rules regarding tDelete, reuse, GetStore
  342. //    All matrices processed during expression evaluation must be subject
  343. //    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
  344. //    If reuse returns TRUE the matrix must be reused.
  345. //    GetMatrix(gm) always calls gm->GetStore()
  346. //    gm->Evaluate obeys rules
  347. //    bm->Evaluate obeys rules for matrices in bm structure
  348.  
  349. void GeneralMatrix::tDelete()
  350. {
  351.    if (tag<0)
  352.    {
  353.       if (tag<-1) { REPORT store=0; delete this; return; }  // borrowed
  354.       else { REPORT return; }
  355.    }
  356.    if (tag==1)
  357.    {
  358.       REPORT  MONITOR_REAL_DELETE("Free   (tDelete)",storage,store)
  359. #ifdef Version21
  360.       if (store) delete [] store;
  361. #else
  362.       if (store) delete [storage] store;
  363. #endif
  364.       store=0; tag=-1; return;
  365.    }
  366.    if (tag==0) { REPORT delete this; return; }
  367.    REPORT tag--; return;
  368. }
  369.  
  370. static void BlockCopy(int n, Real* from, Real* to)
  371. {
  372.    REPORT
  373.    int i = (n >> 3);
  374.    while (i--)
  375.    {
  376.       *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
  377.       *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
  378.    }
  379.    i = n & 7; while (i--) *to++ = *from++;
  380. }
  381.  
  382. Boolean GeneralMatrix::reuse()
  383. {
  384.    if (tag<-1)
  385.    {
  386.       REPORT
  387.       Real* s = new Real [storage]; MatrixErrorNoSpace(s);
  388.       MONITOR_REAL_NEW("Make     (reuse)",storage,s)
  389.       BlockCopy(storage, store, s); store=s; tag=0; return TRUE;
  390.    }
  391.    if (tag<0) { REPORT return FALSE; }
  392.    if (tag<=1)  { REPORT return TRUE; }
  393.    REPORT tag--; return FALSE;
  394. }
  395.  
  396. Real* GeneralMatrix::GetStore()
  397. {
  398.    if (tag<0 || tag>1)
  399.    {
  400.       Real* s;
  401.       if (storage)
  402.       {
  403.          s = new Real [storage]; MatrixErrorNoSpace(s);
  404.          MONITOR_REAL_NEW("Make  (GetStore)",storage,s)
  405.          BlockCopy(storage, store, s);
  406.       }
  407.       else s = 0;
  408.       if (tag>1) { REPORT tag--; }
  409.       else if (tag < -1) { REPORT store=0; delete this; } // borrowed store
  410.       else { REPORT }
  411.       return s;
  412.    }
  413.    Real* s=store; store=0;
  414.    if (tag==0) { REPORT delete this; }
  415.    else { REPORT tag=-1; }
  416.    return s;
  417. }
  418.  
  419. /*
  420. #ifndef __ZTC__
  421. void GeneralMatrix::GetMatrixC(const GeneralMatrix* gmx)
  422. {
  423.    REPORT tag=-1;
  424.    nrows=gmx->nrows; ncols=gmx->ncols; storage=gmx->storage;
  425.    SetParameters(gmx); 
  426.    store = new Real [storage]; MatrixErrorNoSpace(store);
  427.    MONITOR_REAL_NEW("Make (GetMatrix)",storage,store)
  428.    BlockCopy(storage, gmx->store, store);
  429. }
  430. #endif
  431. */
  432.  
  433. void GeneralMatrix::GetMatrix(const GeneralMatrix* gmx)
  434. {
  435.    REPORT  tag=-1; nrows=gmx->Nrows(); ncols=gmx->Ncols();
  436.    storage=gmx->storage; SetParameters(gmx);
  437.    store=((GeneralMatrix*)gmx)->GetStore();
  438. }
  439.  
  440. GeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
  441. // Copy storage of *this to storage of *gmx. Then convert to type mt.
  442. // If mt == 0 just let *gmx point to storage of *this if tag==-1.
  443. {
  444.    if (!mt)
  445.    {
  446.       if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
  447.       else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
  448.    }
  449.    else if (Compare(gmx->Type(),mt))
  450.    { REPORT  gmx->tag = 0; gmx->store = GetStore(); }
  451.    else
  452.    {
  453.       REPORT gmx->tag = -2; gmx->store = store;
  454.       gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
  455.    }
  456.  
  457.    return gmx;
  458. }
  459.  
  460. void GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt)
  461. // Count number of references to this in X.
  462. // If zero delete storage in X;
  463. // otherwise tag X to show when storage can be deleted
  464. // evaluate X and copy to current object
  465. {
  466.    int counter=X.search(this);
  467.    if (counter==0)
  468.    {
  469.       REPORT
  470.       if (store)
  471.       {
  472.          MONITOR_REAL_DELETE("Free (operator=)",storage,store)
  473. #ifdef Version21
  474.          REPORT delete [] store; storage=0;
  475. #else
  476.          REPORT delete [storage] store; storage=0;
  477. #endif
  478.       }
  479.    }
  480.    else { REPORT Release(counter); }
  481.    GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
  482.    if (gmx!=this) { REPORT GetMatrix(gmx); }
  483.    else { REPORT }
  484.    Protect();
  485. }
  486.  
  487. void GeneralMatrix::Inject(const GeneralMatrix& X)
  488. // copy stored values of X; otherwise leave els of *this unchanged
  489. {
  490.    REPORT
  491.    Tracer tr("Inject");
  492.    if (nrows != X.nrows || ncols != X.ncols)
  493.       Throw(IncompatibleDimensionsException());
  494.    MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
  495.    MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
  496.    int i=nrows;
  497.    while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
  498. }  
  499.  
  500. /*************** checking for data loss during conversion *******************/
  501.  
  502. //void GeneralMatrix::CheckConversion(const BaseMatrix& M)
  503. //{
  504. //   if (!(this->Type() >= M.Type()))
  505. //      Throw(ProgramException("Illegal Conversion"));
  506. //}
  507.  
  508. Boolean MatrixConversionCheck::DoCheck;     // will be set to FALSE
  509.  
  510. void MatrixConversionCheck::DataLoss()
  511.    { if (DoCheck) Throw(ProgramException("Illegal Conversion")); }
  512.  
  513. Boolean Compare(const MatrixType& source, MatrixType& destination)
  514. {
  515.    if (!destination) { destination=source; return TRUE; }
  516.    if (destination==source) return TRUE;
  517.    if (MatrixConversionCheck::IsOn() && !(destination>=source))
  518.       Throw(ProgramException("Illegal Conversion", source, destination));
  519.    return FALSE;
  520. }
  521.  
  522. /*************** Make a copy of a matrix on the heap *********************/
  523.  
  524. GeneralMatrix* Matrix::Image() const
  525. {
  526.    REPORT
  527.    GeneralMatrix* gm = new Matrix(*this); MatrixErrorNoSpace(gm);
  528.    return gm;
  529. }
  530.  
  531. GeneralMatrix* SymmetricMatrix::Image() const
  532. {
  533.    REPORT
  534.    GeneralMatrix* gm = new SymmetricMatrix(*this); MatrixErrorNoSpace(gm);
  535.    return gm;
  536. }
  537.  
  538. GeneralMatrix* UpperTriangularMatrix::Image() const
  539. {
  540.    REPORT
  541.    GeneralMatrix* gm = new UpperTriangularMatrix(*this);
  542.    MatrixErrorNoSpace(gm); return gm;
  543. }
  544.  
  545. GeneralMatrix* LowerTriangularMatrix::Image() const
  546. {
  547.    REPORT
  548.    GeneralMatrix* gm = new LowerTriangularMatrix(*this);
  549.    MatrixErrorNoSpace(gm); return gm;
  550. }
  551.  
  552. GeneralMatrix* DiagonalMatrix::Image() const
  553. {
  554.    REPORT
  555.    GeneralMatrix* gm = new DiagonalMatrix(*this); MatrixErrorNoSpace(gm);
  556.    return gm;
  557. }
  558.  
  559. GeneralMatrix* RowVector::Image() const
  560. {
  561.    REPORT
  562.    GeneralMatrix* gm = new RowVector(*this); MatrixErrorNoSpace(gm);
  563.    return gm;
  564. }
  565.  
  566. GeneralMatrix* ColumnVector::Image() const
  567. {
  568.    REPORT
  569.    GeneralMatrix* gm = new ColumnVector(*this); MatrixErrorNoSpace(gm);
  570.    return gm;
  571. }
  572.  
  573. GeneralMatrix* BandMatrix::Image() const
  574. {
  575.    REPORT
  576.    GeneralMatrix* gm = new BandMatrix(*this); MatrixErrorNoSpace(gm);
  577.    return gm;
  578. }
  579.  
  580. GeneralMatrix* UpperBandMatrix::Image() const
  581. {
  582.    REPORT
  583.    GeneralMatrix* gm = new UpperBandMatrix(*this); MatrixErrorNoSpace(gm);
  584.    return gm;
  585. }
  586.  
  587. GeneralMatrix* LowerBandMatrix::Image() const
  588. {
  589.    REPORT
  590.    GeneralMatrix* gm = new LowerBandMatrix(*this); MatrixErrorNoSpace(gm);
  591.    return gm;
  592. }
  593.  
  594. GeneralMatrix* SymmetricBandMatrix::Image() const
  595. {
  596.    REPORT
  597.    GeneralMatrix* gm = new SymmetricBandMatrix(*this); MatrixErrorNoSpace(gm);
  598.    return gm;
  599. }
  600.  
  601. GeneralMatrix* nricMatrix::Image() const
  602. {
  603.    REPORT
  604.    GeneralMatrix* gm = new nricMatrix(*this); MatrixErrorNoSpace(gm);
  605.    return gm;
  606. }
  607.  
  608. GeneralMatrix* GeneralMatrix::Image() const
  609. {
  610.    REPORT
  611.    Throw(InternalException("Cannot apply Image to this matrix type"));
  612.    return 0;
  613. }
  614.  
  615.  
  616. /************************* nricMatrix routines *****************************/
  617.  
  618. void nricMatrix::MakeRowPointer()
  619. {
  620.    row_pointer = new Real* [nrows]; MatrixErrorNoSpace(row_pointer);
  621.    Real* s = Store() - 1; int i = nrows; Real** rp = row_pointer;
  622.    while (i--) { *rp++ = s; s+=ncols; }
  623. }
  624.  
  625. void nricMatrix::DeleteRowPointer()
  626. #ifdef Version21
  627. { if (nrows) delete [] row_pointer; }
  628. #else
  629. { if (nrows) delete [nrows] row_pointer; }
  630. #endif
  631.  
  632. void GeneralMatrix::CheckStore() const
  633. {
  634.    if (!store) 
  635.       Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
  636. }
  637.  
  638.  
  639. /***************************** CleanUp routines *****************************/
  640.  
  641. void GeneralMatrix::CleanUp()
  642. {
  643.    // set matrix dimensions to zero, delete storage
  644.    REPORT
  645.    if (store && storage)
  646.    {
  647.       MONITOR_REAL_DELETE("Free (CleanUp)    ",storage,store)
  648. #ifdef Version21
  649.          REPORT delete [] store;
  650. #else
  651.          REPORT delete [storage] store;
  652. #endif
  653.    }
  654.    store=0; storage=0; nrows=0; ncols=0;
  655. }
  656.  
  657. void nricMatrix::CleanUp()
  658. { DeleteRowPointer(); GeneralMatrix::CleanUp(); }
  659.  
  660. void RowVector::CleanUp()
  661. { GeneralMatrix::CleanUp(); nrows=1; }
  662.  
  663. void ColumnVector::CleanUp()
  664. { GeneralMatrix::CleanUp(); ncols=1; }
  665.  
  666. void CroutMatrix::CleanUp()
  667. { 
  668. #ifdef Version21
  669.    if (nrows) delete [] indx;
  670. #else
  671.    if (nrows) delete [nrows] indx;
  672. #endif
  673.    GeneralMatrix::CleanUp();
  674. }
  675.  
  676. void BandLUMatrix::CleanUp()
  677. { 
  678. #ifdef Version21
  679.    if (nrows) delete [] indx;
  680.    if (storage2) delete [] store2;
  681. #else
  682.    if (nrows) delete [nrows] indx;
  683.    if (storage2) delete [storage2] store2;
  684. #endif
  685.    GeneralMatrix::CleanUp();
  686. }
  687.  
  688.  
  689.